load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

begin
  model = (/"ACCESS-CM2","ACCESS-ESM1-5","BCC-CSM2-MR","CNRM-CM6-1","GFDL-ESM4",\
            "HadGEM3-GC31-LL","IPSL-CM6A-LR","MIROC6","MRI-ESM2-0","NorESM2-LM"/); 10 CMIP6 models

  fread  = addfile("data/pr_eof_1951-2020_JJAS_GPCC.nc","r")
  eof0   = fread->eof_regress(0,:,:)

  fread  = addfile("data/obs/rainfall/GPCC/full_data_monthly_v2022_10.nc","r")
  lat    = fread->lat({25:42})
  lon    = fread->lon({65:105})
  eof    = linint2_Wrap(eof0&lon,eof0&lat,eof0,False,lon,lat,0); regrid to 1°×1°

  fread1 = addfile("data/topo/GMTED2010/GMTED2010_15n030_0125deg.nc","r")
  h      = short2flt(fread1->elevation)*1.
  lat0   = fread1->latitude
  lon0   = fread1->longitude

  h!0    = "lat"
  h!1    = "lon"
  h&lat  = lat0
  h&lon  = lon0
  h      = h/1000.

  hh  = h({25:42},{65:105})
  h  := linint2_Wrap(hh&lon,hh&lat,hh,False,lon,lat,0)
  eof = where(h.ge.2.5,eof,eof@_FillValue); only keep the pattern over HMA

  do i=0,9

    print(model(i))

    dilo  = "data/CMIP6_DA/"+model(i)+"/piControl/mon/pr/"
    fils  = systemfunc("ls "+dilo+"pr_Amon_"+model(i)+"_piControl*.nc")
    fread = addfile(fils,"r")
    pr0   = fread->pr(:,{20:50},{60:110})
    pr0   = pr0*86400*30.; mm/month
    pr1   = linint2_Wrap(pr0&lon,pr0&lat,pr0,False,eof&lon,eof&lat,0); regrid to 1°×1°
    pr    = month_to_season(pr1,"JJA")
    pr    = (pr1(5::12,:,:)+pr1(6::12,:,:)+pr1(7::12,:,:)+pr1(8::12,:,:))/4. ; summer mean
    pr    = where(ismissing(conform(pr,eof,(/1,2/))),pr@_FillValue,pr); only keep the precipitation over HMA

    ts0 = regCoef_n(eof,pr,(/0,1/),(/1,2/)) ;timeseries

    ts0 = dtrend_n(ts0,False,0); remove the overall drift

    n  = dimsizes(ts0)
    nn = n/70 

    ts = new((/nn,70/),"float") 
    do j=0,nn-1 
      ts(j,:) = ts0(j*70:j*70+69);  divided into chunks to have the same length as the observation
    end do 

    foutname = "data/pr_pc1_1951-2020_JJAS_piControl_"+model(i)+".nc"; timeseries for each model
    system("rm -rf "+foutname)
    fout = addfile(foutname,"c")
    fout->ts = ts

    delete([/pr0,pr1,pr,ts0,ts/])
  end do 

end